home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / a_utils / expanded.lha / expanded / src / GeOb.C < prev    next >
C/C++ Source or Header  |  1992-03-19  |  61KB  |  2,228 lines
  1. //
  2. // Linear-Affine-Projective Geometry Package
  3. //
  4. // GeOb.C
  5. //
  6. // $Header$
  7. //
  8. // William J.R. Longabaugh 
  9. // University of Washington
  10. //
  11. // Implementation of the linear-affine-projective geometry
  12. // package described in William J.R. Longabaugh, "An Expanded
  13. // System for Coordinate-Free Geometric Programming", Master's 
  14. // thesis, University of Washington, 1992.
  15. //
  16. // Copyright (c) 1992, William J.R. Longabaugh
  17. //   Copying, use and development for non-commercial purposes permitted.
  18. //   All rights for commercial use reserved.
  19. //   This software is unsupported and without warranty; it is
  20. //   provided "as is".
  21. //
  22. // ***********************************************************************
  23. //
  24. // Notes on efficiency:
  25. //
  26. //   There is a note in the header file Geom.h on efficiency that is
  27. // applicable here.  Unfortunately, the casting of GeObs from one type to
  28. // another is probably very common, and yet the "first cut"
  29. // implementation used here is incredibly inefficient.  Part of the
  30. // problem is that the functions to do the casting (e.g. GeOb
  31. // GeOb::MapTo(GeObType t)) use functions provided to the user (e.g. GeOb
  32. // Map::operator()(GeOb& v)) to do their job.  This means that LOTS of
  33. // redundant checks are being done, and the object is being copied MANY
  34. // times during the course of the casting (thanks to the functional
  35. // semantics of the user functions).  The next cut at improving the
  36. // implementation should provide special functions, available ONLY to the
  37. // implementing classes, that place their results in a designated
  38. // location and avoid type checking/casting their arguments, e.g.
  39. //       GeOb Map::Apply(GeOb& v, GeOb* out);
  40. // This approach would help a lot at making things more efficient.
  41. //
  42. //   Lots of implementing functions use the approach to call CanMap()
  43. // first to check if a map will succeed, and then calling MapTo().  They
  44. // do this to make error messages to the user more meaningful, since
  45. // errors will be reported by the function the user calls, and not a
  46. // function lower down.  However, this is very inefficient, since it
  47. // means the GeOb is being checked if it is OK more than once.  A better
  48. // approach would have GeObs providing the function
  49. //       Boolean GeOb::MapTo(GeObType t, GeOb* out);
  50. // (ONLY to the implementing classes) that would avoid object copying
  51. // and return FALSE if the mapping could not be done.
  52. //
  53. // ***********************************************************************
  54.  
  55. #include <math.h>
  56. #include "Lap.h"
  57.  
  58. static Scalar randval(void);
  59.  
  60. // ***********************************************************************
  61. //
  62. // This function returns a NON-ZERO integer in the range -10.0 to +10.0,
  63. // cast to a Scalar:
  64. //
  65.  
  66. #define RAND_RANGE 20
  67.  
  68. static Scalar randval()
  69. {
  70. #ifdef c_plusplus
  71.   long random();
  72. #endif
  73.  
  74.   Scalar retval;
  75.  
  76.   while (TRUE) { 
  77.     retval = (Scalar)((random() % RAND_RANGE) - (RAND_RANGE / 2));
  78.     if (retval != 0.0) break;
  79.   }
  80.   return retval;
  81. }
  82.  
  83. // ***********************************************************************
  84. // ***********************************************************************
  85. //
  86. // GeOb Class
  87. //
  88. // ***********************************************************************
  89. // ***********************************************************************
  90. //
  91. // Private/protected member functions
  92. //
  93. // ***********************************************************************
  94. //
  95. // Used by various classes to create GeObs with given matrix
  96. // representation:
  97. //
  98.  
  99. GeOb::GeOb(GeObType g, Space& ins, Matrix& inm) : s(ins), m(inm)
  100. {
  101.   type = ANY_GEOB;
  102.   holds = g;
  103. }
  104.  
  105. // ***********************************************************************
  106. //
  107. // Dumps data for GeObs:
  108. //
  109.  
  110. void GeOb::data_out(ostream &c, int indent, char* name)
  111. {
  112.   char *ibloc = new char[indent + 1];
  113.   for (int i = 0; i < indent; i++) {
  114.     *(ibloc + i) = ' ';
  115.   }
  116.   *(ibloc + indent) = '\0';
  117.  
  118.   c << ibloc << ast;
  119.   c << ibloc << name;
  120.   c << ibloc << "Type of GeOb: ";
  121.   GeObTypeOut(c, type);
  122.   c << "\n";
  123.   c << ibloc << "Currently holds: ";
  124.   GeObTypeOut(c, holds);
  125.   c << "\n";
  126.   if (holds != NULL_GEOB) {
  127.     c << ibloc << "Contained in space:\n";
  128.     s.debug_out(c, indent + ERR_IND);
  129.     c << ibloc << "Matrix representation of the object:\n";
  130.     m.debug_out(c, indent + ERR_IND);
  131.   }
  132.   c << ibloc << ast;
  133.  
  134.   delete ibloc;
  135.   return;
  136. }
  137.  
  138. // ***********************************************************************
  139. //
  140. // Public member functions
  141. //
  142. // ***********************************************************************
  143. //
  144. // Need to do memberwise initialization, since Matrix class members need
  145. // to do memberwise initialization:
  146. //
  147.  
  148. GeOb::GeOb(GeOb& g) : s(g.s), m(g.m)
  149. {
  150.   type = ANY_GEOB;
  151.   holds = g.holds;
  152. }
  153.  
  154. // ***********************************************************************
  155. //
  156. // Need to do memberwise assignment, since Matrix class members need
  157. // to do memberwise assignment:
  158. //
  159.  
  160. GeOb& GeOb::operator=(GeOb& g)
  161. {
  162. // If there is a mismatch between what the arg. holds and the type
  163. // of this Geob, and this geob is not an ANY_GEOB, we need to cast.
  164. // This checking needs to be done because of apparent assignment
  165. // inheritance in g++ 1.37:
  166.  
  167.   if ((type != ANY_GEOB) && (type != g.holds) && (g.holds != NULL_GEOB)) {
  168.     *this = g.MapTo(type);
  169.   } else {
  170.     holds = g.holds;
  171.     s = g.s;
  172.     m = g.m;
  173.   }
  174.   return (*this);
  175. }
  176.  
  177. // ***********************************************************************
  178. //
  179. // Output for debugging:
  180. //
  181.  
  182. void GeOb::debug_out(ostream &c, int indent)
  183. {
  184.   static char* name = "GeOb Object\n";
  185.   this->data_out(c, indent, name);
  186.   return;
  187. }
  188.  
  189. // ***********************************************************************
  190. //
  191. // If the multimap has zero arguments, we can cast it into a GeOb in
  192. // the range space.  If it is a first-order map into the reals, it can
  193. // be cast into a vector in the dual space of the domain.
  194. //
  195.  
  196. GeOb::GeOb(MultiMap& mp)
  197. {
  198.   Boolean success = FALSE;
  199.  
  200.   type = ANY_GEOB;
  201.  
  202. // Find out if it has zero arguments.  If it does, do the cast.
  203. // Get the matrix rep. of the object.  It is possible, if the type
  204. // is AFF_VECTOR, that the representation needs to be de-hoisted:
  205.  
  206.   if (mp.FullyEvaluated()) {
  207.     m = mp.FullEvalRep();
  208.     s = mp.RangeSpace();
  209.     holds = s.NativeType();
  210.     if ((holds == AFF_VECTOR) && (m.Columns() != s.MatrixSize())) {
  211.       m *= Transpose(s.GetSpace(AFFINE).HoistTangent());
  212.     }
  213.     success = TRUE;
  214.   } else {
  215.  
  216. //  Otherwise, check if it is a linear functional:
  217.  
  218.     Space range(mp.RangeSpace());
  219.     Space domain(mp.DomainSpace());
  220.  
  221.     if ((mp.Holds() == MULTI_LINEAR) && (range == Reals) &&
  222.     (domain.CPSpaceSize() == 1)) {
  223.       s = domain.Dual();
  224.       m = Matrix(1, s.MatrixSize());
  225.       holds = s.NativeType();
  226.       for (int i = 0; i < s.MatrixSize(); i++) {
  227.         m[0][i] = (mp.GetStdImage(i))[0][0];
  228.       }
  229.       success = TRUE;
  230.     }
  231.   }
  232.  
  233. // If it is neither, we have an error:
  234.  
  235.   if (!success) {
  236.     errh.ErrorExit("GeOb::GeOb(MultiMap&)",
  237.                    "MultiMap cannot be cast to a GeOb", mp);
  238.   }
  239. }
  240.  
  241. // ***********************************************************************
  242. //
  243. // Casts maps that are unrestricted linear functionals into dual vectors:
  244. //
  245.  
  246. GeOb::GeOb(Map& mp)
  247. {
  248. // Check the range and domain of the map.  If the domain is not a
  249. // subset, and the range is the Reals, the map can be converted
  250. // into a Vector in the dual space of the domain:
  251.  
  252.   SubSet Range(mp.Range());
  253.   SubSet Domain(mp.Domain());
  254.   Space  embeds(Range.EmbeddingSpace());
  255.  
  256.   type = ANY_GEOB;
  257.  
  258.   if ((mp.Holds() == LIN_MAP) && (embeds == Reals) &&
  259.       Range.IsFullSpace() && Domain.IsFullSpace()) {
  260.     s = (Domain.EmbeddingSpace()).Dual();
  261.     m = Transpose(mp.MatrixRep());
  262.     holds = s.NativeType();
  263.   } else {
  264.     errh.ErrorExit("GeOb::GeOb(Map&)",
  265.                    "Map is not an unrestricted linear functional", mp);
  266.   }  
  267. }
  268.  
  269. // ***********************************************************************
  270. //
  271. // Casts Scalars into vectors in the predefined space "Reals", using
  272. // the standard basis, i.e. "1":
  273. //
  274.  
  275. GeOb::GeOb(Scalar v) : s(Reals), m(IdentityMatrix(1) * v)
  276. {
  277.   type = ANY_GEOB;
  278.   holds = VECTOR;
  279. }
  280.  
  281. // ***********************************************************************
  282. //
  283. // Says if the GeOb is (or can be mapped to) a zero vector:
  284. //
  285.  
  286. Boolean GeOb::IsZeroVector(void)
  287. {
  288.   Matrix mat = m;
  289.  
  290.   if ((holds != VECTOR) && (holds != AFF_VECTOR)) {
  291.     GeOb temp = this->MapTo(VECTOR);
  292.     mat = temp.m;
  293.   }
  294.  
  295.   for (int i = 0; i < s.Dim(); i++) {
  296.     if (fabs(mat[0][i]) > EPSILON) {
  297.       return (FALSE);
  298.     }
  299.   }
  300.   return (TRUE);
  301. }
  302.  
  303. // ***********************************************************************
  304. //
  305. // This returns the nth element of a vector or point tuple (a vector
  306. // or point in a cartesian product space).
  307. //
  308.  
  309. GeOb GeOb::operator[](int n)
  310. {
  311.   GeOb temp(*this);
  312.  
  313. // By storing points using a standard frame, point
  314. // extraction is pretty straightforward.  If we were to
  315. // use a simplex, we would need to sum some coordinates to get
  316. // the coordinates to return.
  317.  
  318. // This only works if this GeOb is a vector, aff_vector, or point.
  319. // If it is not, map it into the linearization space.  No other attempt
  320. // will be made to "find" a space that is a CP space.  (The user may need to
  321. // explicitly apply standard maps).
  322.  
  323.   if ((holds != VECTOR) && (holds != AFF_VECTOR) && (holds != AFF_POINT)) {
  324.     temp = this->MapTo(VECTOR);
  325.   }
  326.  
  327. // Check that the integer is in the correct range:
  328.  
  329.   Space tempsp(temp.SpaceOf());
  330.   if ((n >= tempsp.CPSpaceSize()) || (n < 0)) {
  331.     errh.ErrorExit("GeOb GeOb::operator[](int)",
  332.                    "Specified n is invalid",
  333.                    ErrVal("Value of n = ", n), temp);
  334.   }
  335.  
  336. // If the mapped object is an atomic object, just return it:
  337.  
  338.   if (tempsp.CPSpaceSize() == 1) {
  339.     return (temp);
  340.   }
  341.  
  342. // Get the nth space in the tuple and its coordinate starting slot,
  343. // and figure out the size of the matrix to build:
  344.  
  345.   Space retspace = tempsp[n];
  346.   int slot = tempsp.StartSlot(n);
  347.  
  348. // Build the matrix.  If this is a point, we need to tack on the
  349. // trailing 1:
  350.  
  351.   Matrix retmat(1, retspace.MatrixSize());
  352.   for (int i = 0; i < retspace.Dim(); i++) {
  353.     retmat[0][i] = (temp.m)[0][slot++];
  354.   }
  355.   if (retspace.Holds() == AFF_SPACE) {
  356.     retmat[0][retspace.Dim()] = 1.0;
  357.   }
  358.  
  359.   GeOb retval(retspace.NativeType(), retspace, retmat);
  360.   return (retval);
  361. }
  362.  
  363. // ***********************************************************************
  364. //
  365. // Like the above function, but it extracts all the elements from
  366. // the tuple and returns them in a list:
  367.  
  368. GeObList GeOb::TupleElements(void)
  369. {
  370.   GeOb temp(*this);
  371.  
  372. // Same comments apply here as with above function.
  373.  
  374.   if ((holds != VECTOR) && (holds != AFF_VECTOR) && (holds != AFF_POINT)) {
  375.     temp = this->MapTo(VECTOR);
  376.   }
  377.  
  378. // Loop through all the elements, creating a new GeOb for each element,
  379. // and adding it to the list of GeObs that is returned.
  380. // (But if the mapped object is an atomic object, just return it):
  381.  
  382.   Space tempsp(temp.SpaceOf());
  383.   int n = tempsp.CPSpaceSize();
  384.   GeObList retlist(n);
  385.  
  386.   if (tempsp.CPSpaceSize() == 1) {
  387.     retlist[0] = temp;
  388.     return (retlist);
  389.   }
  390.  
  391.   for (int j = 0; j < n; j++) {
  392.     Space retspace = tempsp[j];
  393.     int slot = tempsp.StartSlot(j);
  394.  
  395. //  Build the matrix.  If this is a point, we need to tack on the
  396. //  trailing 1:
  397.  
  398.     Matrix retmat(1, retspace.MatrixSize());
  399.     for (int i = 0; i < retspace.Dim(); i++) {
  400.       retmat[0][i] = (temp.m)[0][slot++];
  401.     }
  402.     if (retspace.Holds() == AFF_SPACE) {
  403.       retmat[0][retspace.Dim()] = 1.0;
  404.     }
  405.  
  406.     retlist[j] = GeOb(retspace.NativeType(), retspace, retmat);
  407.   }
  408.   return (retlist);
  409. }
  410.  
  411. // ***********************************************************************
  412. //
  413. // IMPORTANT NOTE ON ALGEBRAIC AND APPLICATION OPERATIONS:
  414. //
  415. // In the following operations, it is considered to make sense to
  416. // use these operators on vector equivalence classes (ec) (and affine
  417. // vector equivalence classes).  Since casting ec to vectors
  418. // or affine vectors returns random elements, the result of the
  419. // operations will also be typically random.  While it is not
  420. // recommended to use the operations on ec, the operations are
  421. // included for completeness.
  422. //
  423. // ***********************************************************************
  424. //
  425. // FRIEND FUNCTION
  426. // If the object is a vector, negate it.  Otherwise, map it into
  427. // the linearization space and negate it there (though try to keep
  428. // operations in tangent spaces localized there).
  429. //
  430.  
  431. GeOb operator-(GeOb& th)
  432. {
  433.   if ((th.holds == VECTOR) || (th.holds == AFF_VECTOR)) {
  434.     GeOb retval(th.holds, th.s, -(th.m));
  435.     return (retval);
  436.   } else if (th.holds == AFF_VEC_EC) {
  437.     GeOb retval = th.MapTo(AFF_VECTOR);
  438.     retval.m = -(retval.m);
  439.     return (retval);
  440.   } else {
  441.     GeOb retval = th.MapTo(VECTOR);
  442.     retval.m = -(retval.m);
  443.     return (retval);
  444.   }
  445. }
  446.  
  447. // ***********************************************************************
  448. //
  449. // FRIEND FUNCTION
  450. // If both objects are vectors in the same space, add them.  Try to
  451. // keep operations in tangent spaces localized there.  Otherwise,
  452. // map objects to their Linearization spaces and add them there if the
  453. // spaces match.
  454. //
  455.  
  456. GeOb operator+(GeOb& th, GeOb& g)
  457. {
  458.   static char* sig = "GeOb operator+(GeOb&, GeOb&)";
  459.  
  460.   if (((th.holds == VECTOR) && (g.holds == VECTOR)) ||
  461.       ((th.holds == AFF_VECTOR) && (g.holds == AFF_VECTOR))) {
  462.     if (th.s == g.s) {
  463.       GeOb retval(th.holds, th.s, th.m + g.m);
  464.       return (retval);
  465.     } else {
  466.       errh.ErrorExit(sig, "Spaces of operands do not match", th, g);
  467.     }
  468.   } else if (((th.holds == AFF_VECTOR) || (th.holds == AFF_VEC_EC)) &&
  469.              ((g.holds == AFF_VEC_EC) || (g.holds == AFF_VECTOR))) {
  470.     GeOb retval = th.MapTo(AFF_VECTOR);
  471.     GeOb other = g.MapTo(AFF_VECTOR);
  472.     if (retval.s == other.s) {
  473.       retval.m = retval.m + other.m;
  474.       return (retval);
  475.     } else {
  476.       errh.ErrorExit(sig, "Spaces of operands do not match", th, g);
  477.     }
  478.   } else {
  479.     GeOb retval = th.MapTo(VECTOR);
  480.     GeOb other = g.MapTo(VECTOR);
  481.     if (retval.s == other.s) {
  482.       retval.m = retval.m + other.m;
  483.       return (retval);
  484.     } else {
  485.       errh.ErrorExit(sig, "Spaces of mapped operands do not match",
  486.                      th, g, retval, other);
  487.     }
  488.   }
  489. }
  490.  
  491. // ***********************************************************************
  492. //
  493. // FRIEND FUNCTION
  494. // Handle subtraction just like addition:
  495. //
  496.  
  497. GeOb operator-(GeOb& th, GeOb& g)
  498. {
  499.   static char* sig = "GeOb operator-(GeOb&, GeOb&)";
  500.  
  501.   if (((th.holds == VECTOR) && (g.holds == VECTOR)) ||
  502.       ((th.holds == AFF_VECTOR) && (g.holds == AFF_VECTOR))) {
  503.     if (th.s == g.s) {
  504.       GeOb retval(th.holds, th.s, th.m - g.m);
  505.       return (retval);
  506.     } else {
  507.       errh.ErrorExit(sig, "Spaces of operands do not match", th, g);
  508.     }
  509.   } else if (((th.holds == AFF_VECTOR) || (th.holds == AFF_VEC_EC)) &&
  510.              ((g.holds == AFF_VEC_EC) || (g.holds == AFF_VECTOR))) {
  511.     GeOb retval = th.MapTo(AFF_VECTOR);
  512.     GeOb other = g.MapTo(AFF_VECTOR);
  513.     if (retval.s == other.s) {
  514.       retval.m = retval.m - other.m;
  515.       return (retval);
  516.     } else {
  517.       errh.ErrorExit(sig, "Spaces of operands do not match", th, g);
  518.     }
  519.   } else {
  520.     GeOb retval = th.MapTo(VECTOR);
  521.     GeOb other = g.MapTo(VECTOR);
  522.     if (retval.s == other.s) {
  523.       retval.m = retval.m - other.m;
  524.       return (retval);
  525.     } else {
  526.       errh.ErrorExit(sig, "Spaces of mapped operands do not match",
  527.                      th, g, retval, other);
  528.     }
  529.   }
  530. }
  531.  
  532. // ***********************************************************************
  533. //
  534. // FRIEND FUNCTION
  535. // This operation only make sense if one of the operands is a Vector
  536. // in the special 1-dim space "Reals", in which case this is scalar
  537. // multiplication.  While it would be very nice to be able to define
  538. // functions like "GeOb operator*(Scalar coeff, GeOb& g)" for efficient
  539. // scalar multiplication, compiler ambiguity prevents this from being
  540. // possible:
  541. //
  542.  
  543. GeOb operator*(GeOb& th, GeOb& g)
  544. {
  545.   GeOb temp;
  546.   Scalar coeff;
  547.  
  548.   if (th.s == Reals) {
  549.     coeff = th.ToScalar();
  550.     temp = g;
  551.   } else if (g.s == Reals) {
  552.     coeff = g.ToScalar();
  553.     temp = th;
  554.   } else {
  555.     errh.ErrorExit("GeOb operator*(GeOb&, GeOb&)",
  556.                    "Neither operand was a vector in space Reals", th, g);
  557.   }
  558.  
  559. // Try to keep operations local to tangent vector spaces:
  560.  
  561.   if (temp.holds == AFF_VEC_EC) {
  562.     temp = temp.MapTo(AFF_VECTOR);
  563.   } else if ((temp.holds != VECTOR) && (temp.holds != AFF_VECTOR)) {
  564.     temp = temp.MapTo(VECTOR);
  565.   }
  566.  
  567.   temp.m = temp.m * coeff;
  568.   return (temp);
  569. }
  570.  
  571. // ***********************************************************************
  572. //
  573. // FRIEND FUNCTION
  574. // Same as multiplication, using a reciprocal.  The second argument
  575. // must be the vector from the Reals, and it must be non-zero.  Again
  576. // having a function "GeOb operator/(GeOb& g, Scalar coeff)" would
  577. // be great, but must be avoided to prevent compiler ambiguity:
  578. //
  579.  
  580. GeOb operator/(GeOb& th, GeOb& g)
  581. {
  582.   static char* sig = "GeOb operator/(GeOb&, GeOb&)";
  583.   Scalar coeff;
  584.  
  585.   if (g.s == Reals) {
  586.     coeff = g.ToScalar();
  587.   } else {
  588.     errh.ErrorExit(sig, "Second operand was not a vector in space Reals",
  589.                    th, g);
  590.   }
  591.  
  592.   if (fabs(coeff) < EPSILON) {
  593.     errh.ErrorExit(sig, "Attempt to divide by zero",
  594.                    ErrVal("Value of denominator = ", coeff), th);
  595.   }
  596.  
  597. // Try to keep operations local to tangent vector spaces:
  598.  
  599.   if ((th.holds == VECTOR) || (th.holds == AFF_VECTOR)) {
  600.     GeOb retval(th.holds, th.s, th.m / coeff);
  601.     return (retval);
  602.   } else if (th.holds == AFF_VEC_EC) {
  603.     GeOb retval = th.MapTo(AFF_VECTOR);
  604.     retval.m = retval.m / coeff;
  605.     return (retval);
  606.   } else {
  607.     GeOb retval = th.MapTo(VECTOR);
  608.     retval.m = retval.m / coeff;
  609.     return (retval);
  610.   }
  611. }
  612.  
  613. // ***********************************************************************
  614. //
  615. // Apply this object to the argument, which should be in the dual space.
  616. // Try to keep operations local to tangent vector spaces if possible.
  617. //
  618.  
  619. Scalar GeOb::Apply(GeOb& a)
  620. {
  621.   static char* sig = "Scalar GeOb::Apply(GeOb&)";
  622.  
  623. // One of the objects is going to be in a dual space, and cannot
  624. // be mapped (except maybe from an equivalence class).  The other
  625. // object should be mapped to the native type of the dual space.
  626.  
  627.   GeObType firsttype;
  628.   GeObType othertype;
  629.   SRel thisrel = s.ThisSpaceIs();
  630.   SRel otherrel = a.SpaceOf().ThisSpaceIs();
  631.   if ((thisrel == TANG_DUAL) || (thisrel == LIN_DUAL)) {
  632.     firsttype = s.NativeType();
  633.     othertype = s.Dual().NativeType();
  634.   } else if ((otherrel == TANG_DUAL) || (otherrel == LIN_DUAL)) {
  635.     firsttype = a.SpaceOf().Dual().NativeType();
  636.     othertype = a.SpaceOf().NativeType();
  637.   } else {
  638.     errh.ErrorExit(sig, "One item must be in a dual space", *this, a);
  639.   }
  640.   GeOb first = this->MapTo(firsttype);
  641.   GeOb other = a.MapTo(othertype);
  642.   if (first.SpaceOf() == (other.SpaceOf()).Dual()) {
  643.     return ((first.m * Transpose(other.m))[0][0]);
  644.   } else {
  645.     errh.ErrorExit(sig, "Spaces of mapped operands are not duals",
  646.                    *this, a, first, other);
  647.   }
  648. }
  649.  
  650. // ***********************************************************************
  651. //
  652. // If the GeOb is a vector in a space that has a distinguished inner
  653. // product (i.e. it is a Euclidean Vector Space), there is a
  654. // mapping between vectors and dual vectors: v -> < ,v>.
  655. //
  656.  
  657. GeOb GeOb::Dual(void)
  658. {
  659.   static char* sig = "Vector GeOb::Dual(void)";
  660.   GeObList hold(2);
  661.  
  662. // This implementation uses the robust way to find the dual: partially
  663. // evaluate the inner product MLM using this vector, then cast the
  664. // resulting linear functional into the dual vector.  While taking the
  665. // transpose is MUCH faster, future changes to the system that permit
  666. // arbitrary distinguished products would break the transpose route.
  667.  
  668.   if ((holds == VECTOR) || (holds == AFF_VECTOR)) {
  669.     if (s.IsEuclidean()) {
  670.       hold[0] = *this;
  671.       return ((s.InnerProduct())(s, hold));
  672.     } else {
  673.       errh.ErrorExit(sig, "Vector in not in a Euclidean space", *this);
  674.     }
  675.   }
  676. // If this GeOb is an affine vector ec, cast it to an affine vector.
  677. // Otherwise cast it to a vector and go:
  678.  
  679.   GeObType maptype = ((holds == AFF_VEC_EC) ? AFF_VECTOR : VECTOR);
  680.   GeOb mapped = this->MapTo(maptype);
  681.   if ((mapped.SpaceOf()).IsEuclidean()) {
  682.     hold[0] = mapped;
  683.     return ((mapped.SpaceOf().InnerProduct())(mapped.SpaceOf(), hold));
  684.   } else {
  685.     errh.ErrorExit(sig, "Mapped object not in a Euclidean space",
  686.            *this, mapped);
  687.   }
  688. }
  689.  
  690. // ***********************************************************************
  691. //
  692. // Returns TRUE iff this object can be cast to the specified object.
  693. // Code is so specific to each type of GeOb that it was decided to
  694. // cast this GeOb down and let functions tailored for each type do
  695. // the work.  While this localizes the implementation, the use of
  696. // downcasting is TERRIBLE for efficiency, since it means this
  697. // GeOb is being recopied prior to the operation being performed:
  698. //
  699.  
  700. Boolean GeOb::CanMapTo(GeObType t)
  701. {
  702.   Boolean rv;
  703.  
  704. // Quick kill:
  705.  
  706.   if (t == holds) return (TRUE);
  707.  
  708. // Need to cover all cases.  Need to ask if this object is in
  709. // the domain of the standard mapping to the given type. 
  710.  
  711.   switch (holds) {
  712.     case VECTOR:
  713.       rv = Vector(*this).CanMapTo(t);
  714.       break;
  715.     case AFF_POINT:
  716.       rv = APoint(*this).CanMapTo(t);
  717.       break;
  718.     case AFF_VECTOR:
  719.       rv = AVector(*this).CanMapTo(t);
  720.       break;
  721.     case VEC_EC:
  722.       rv = VectorEC(*this).CanMapTo(t);
  723.       break;
  724.     case AFF_VEC_EC:
  725.       rv = AVectorEC(*this).CanMapTo(t);
  726.       break;
  727.     case PROJ_POINT:
  728.       rv = PPoint(*this).CanMapTo(t);
  729.       break;
  730.     case NULL_GEOB:
  731.     case ANY_GEOB:
  732.     default:
  733.       errh.ErrorExit("Boolean GeOb::CanMapTo(GeObType)", "Illegal GeOb type",
  734.                  ErrType("Type found =", holds, GEOBTYPES), *this);
  735.       break;
  736.   }
  737.   return (rv);
  738. }
  739.  
  740. // ***********************************************************************
  741. //
  742. // This is the way to get the GeOb that results from standard mappings.
  743. // Code is so specific to each type of GeOb that it was decided to
  744. // cast this GeOb down and let functions tailored for each type do
  745. // the work.  Unfortunately, while this approach localizes the
  746. // implementation, the use of downcasting is TERRIBLE for efficiency,
  747. // since it involves this GeOb being recopied prior to the operation 
  748. // being performed (in fact, it means MapTo() is being called AGAIN,
  749. // though it only executes the "quick kill" base case).  This is
  750. // particulary harmful for system efficiency given the frequency of
  751. // mapping in many common operations (map application, algebraic
  752. // operations).  Another approach may be more desirable.
  753. //
  754.  
  755. GeOb GeOb::MapTo(GeObType t)
  756. {
  757.   GeOb rv;
  758.  
  759. // Quick kill:
  760.  
  761.   if (t == holds) return (*this);
  762.  
  763. // Need to cover all cases.
  764.  
  765.   switch (holds) {
  766.     case VECTOR:
  767.       rv = Vector(*this).MapTo(t);
  768.       break;
  769.     case AFF_POINT:
  770.       rv = APoint(*this).MapTo(t);
  771.       break;
  772.     case AFF_VECTOR:
  773.       rv = AVector(*this).MapTo(t);
  774.       break;
  775.     case VEC_EC:
  776.       rv = VectorEC(*this).MapTo(t);
  777.       break;
  778.     case AFF_VEC_EC:
  779.       rv = AVectorEC(*this).MapTo(t);
  780.       break;
  781.     case PROJ_POINT:
  782.       rv = PPoint(*this).MapTo(t);
  783.       break;
  784.     case NULL_GEOB:
  785.     case ANY_GEOB:
  786.     default:
  787.       errh.ErrorExit("GeOb GeOb::MapTo(GeObType)", "Illegal GeOb type",
  788.              ErrType("Type found =", holds, GEOBTYPES), *this);
  789.       break;
  790.   }
  791.   return (rv);
  792. }
  793.  
  794. // ***********************************************************************
  795. //
  796. // It would be nice to be able to SET tuple elements using the []
  797. // operator, but this doesn't work since a tuple is not implemented
  798. // as a set of GeObs, and so we cannot return a reference to one of them.
  799. //
  800.  
  801. GeOb GeOb::SetTupleElement(int n, GeOb& g)
  802. {
  803.   static char* sig = "GeOb GeOb::SetTupleElement(int, GeOb&)";
  804.   GeOb temp(*this);
  805.  
  806. // First, this only works if this GeOb is a vector, aff_vector, or point.
  807. // If it is not, map it into the linearization space.  No other attempt
  808. // will be made to "find" a space that is a CP space.  (The user may need to
  809. // explicitly apply standard maps).
  810.  
  811.   if ((holds != VECTOR) && (holds != AFF_VECTOR) && (holds != AFF_POINT)) {
  812.     temp = this->MapTo(VECTOR);
  813.   }
  814.  
  815. // Check that the integer is in the correct range:
  816.  
  817.   Space tempsp(temp.SpaceOf());
  818.   if ((n >= tempsp.CPSpaceSize()) || (n < 0)) {
  819.     errh.ErrorExit(sig, "Specified n is invalid",
  820.                    ErrVal("Value of n = ", n), temp);
  821.   }
  822.  
  823. // If the mapped object is in an atomic space, setting the element
  824. // reduces to returning the mapped argument GeOb, as long as the
  825. // spaces match:
  826.  
  827.   if (tempsp.CPSpaceSize() == 1) {
  828.     GeOb retval = g.MapTo(tempsp.NativeType());
  829.     if (tempsp == retval.SpaceOf()) {
  830.       return (retval);
  831.     } else {
  832.       errh.ErrorExit(sig, "Mapped GeOb is not in required space",
  833.                      temp, tempsp, g, retval);
  834.     }
  835.   }
  836.  
  837. // Get the nth space in the tuple and its coordinate starting slot.
  838. // Map the object to the necessary type:
  839.  
  840.   Space slotspace = tempsp[n];
  841.   int slot = tempsp.StartSlot(n);
  842.   GeOb insert = g.MapTo(slotspace.NativeType());
  843.  
  844. // Check that the spaces match, and if they do, replace the
  845. // coordinates in the required slots and return:
  846.  
  847.   if (slotspace == insert.SpaceOf()) {
  848.     Matrix retmat(temp.m); 
  849.     for (int i = 0; i < slotspace.Dim(); i++) {
  850.       retmat[0][slot++] = (insert.m)[0][i];
  851.     }
  852.     GeOb retval = GeOb(temp.Holds(), tempsp, retmat);
  853.     return (retval);
  854.   } else {
  855.     errh.ErrorExit(sig, "Mapped GeOb is not in required space",
  856.                    temp, slotspace, g, insert);
  857.   }
  858. }
  859.  
  860. // ***********************************************************************
  861. //
  862. // FRIEND FUNCTION
  863. // Correspondence of GeObList is INFINITY->A, 0->B, 1->C.
  864. // Given CR = (A,B;C,D) and list (A,B,C), this returns D:
  865. //
  866.  
  867. PPoint CrossRatio(AugScalar v, GeObList& g)
  868. {
  869.   static char* sig = "GeOb GeOb::CrossRatio(AugScalar, GeObList&)";
  870.   int i;
  871.  
  872. // The GeObList must consist of 3 objects that can be cast to
  873. // co-linear points in the projective completion.
  874.  
  875.   if (g.Length() != 3) {
  876.     errh.ErrorExit(sig, "List must contain 3 items", g);
  877.   }
  878.  
  879. // Cast the first GeOb and glean inportant info.
  880.  
  881.   GeOb temp = g[0].MapTo(PROJ_POINT);
  882.   Space retsp = temp.SpaceOf();
  883.   int numcol = retsp.MatrixSize();
  884.   Matrix holdmat(2, numcol);
  885.   Matrix unitmat(1, numcol);
  886.   holdmat[0] = (temp.m)[0]; 
  887.  
  888. // Go through the rest of the list, casting and storing the
  889. // matrix representation:
  890.  
  891.   for (i = 1; i < 3; i++) {
  892.     temp = g[i].MapTo(PROJ_POINT);
  893.     if (temp.SpaceOf() != retsp) {
  894.       errh.ErrorExit(sig, "All items do not map to same space", g);
  895.     }
  896.     if (i == 1) {
  897.       holdmat[1] = (temp.m)[0]; 
  898.     } else {
  899.       unitmat[0] = (temp.m)[0]; 
  900.     }
  901.   }
  902.  
  903. // We need to make sure that the points are distinct and that they
  904. // are all co-linear.  Then we need to scale the rows
  905. // of the matrix so they sum to the unit point matrix.
  906.  
  907. // Build an orthogonal basis using the first two vectors:
  908.  
  909.   Matrix orthbas(2, numcol);
  910.   Scalar mag = sqrt((holdmat[0] * Transpose(holdmat[0]))[0][0]);
  911.   orthbas[0] = (holdmat[0] / mag)[0];
  912.   Scalar fac = (holdmat[1] * Transpose(orthbas[0]))[0][0];
  913.   orthbas[1] = (holdmat[1] - (fac * orthbas[0]))[0];
  914.   mag = sqrt((orthbas[1] * Transpose(orthbas[1]))[0][0]);
  915.   if (fabs(mag) < EPSILON) {
  916.     errh.ErrorExit(sig, "First two points are not distinct", g);
  917.   }
  918.   orthbas[1] = (orthbas[1] / mag)[0];
  919.   
  920. // Project the third point into the plane, and make sure this projection
  921. // is essentially superfluous:
  922.  
  923.   Matrix newunit = unitmat * Transpose(orthbas);
  924.  
  925.   Matrix test = (newunit * orthbas) - unitmat;
  926.   for (i = 0; i < numcol; i++) {
  927.     if (fabs(test[0][i]) > EPSILON) { 
  928.       errh.ErrorExit(sig, "Points are not collinear", g);
  929.     }
  930.   }
  931.  
  932. // Find how to scale the rows of the mapping matrix so they sum to 
  933. // the unit point:
  934.   
  935.   Matrix coeff = newunit * Inverse(holdmat * Transpose(orthbas));
  936.   if ((fabs(coeff[0][0]) < EPSILON) || (fabs(coeff[0][1]) < EPSILON)) { 
  937.     errh.ErrorExit(sig, "Third point is not distinct", g);
  938.   }
  939.   holdmat[0] = (holdmat[0] * coeff[0][0])[0];
  940.   holdmat[1] = (holdmat[1] * coeff[0][1])[0];
  941.  
  942. // That's it!  Holdmat now holds the map from extended reals to
  943. // all points on the line, so map and return:
  944.  
  945.   GeOb retval(PROJ_POINT, retsp, v.GetMatrix() * holdmat);
  946.   return PPoint(retval);
  947. }
  948.  
  949. // ***********************************************************************
  950. //
  951. // FRIEND FUNCTION
  952. // Given a list of four collinear points (A,B,C,D) where at least
  953. // three are distinct, this routine calculates the cross ratio
  954. // (A,B;C,D).
  955. //
  956.  
  957. AugScalar CrossRatioCalc(GeObList& g)
  958. {
  959.   static char* sig = "AugScalar GeOb::CrossRatioCalc(GeObList&)";
  960.   int i;
  961.  
  962. // The GeObList must consist of 4 objects that can be cast to
  963. // co-linear points in the projective completion:
  964.  
  965.   if (g.Length() != 4) {
  966.     errh.ErrorExit(sig, "List must contain 4 items", g);
  967.   }
  968.  
  969. // Cast the first GeOb and glean inportant info.
  970.  
  971.   GeOb temp = g[0].MapTo(PROJ_POINT);
  972.   Space chksp = temp.SpaceOf();
  973.   int numcol = chksp.MatrixSize();
  974.   Matrix holdmat(4, numcol);
  975.   holdmat[0] = (temp.m)[0]; 
  976.  
  977. // Go through the rest of the list, casting and storing the
  978. // matrix representation:
  979.  
  980.   for (i = 1; i < 4; i++) {
  981.     temp = g[i].MapTo(PROJ_POINT);
  982.     if (temp.SpaceOf() != chksp) {
  983.       errh.ErrorExit(sig, "All items do not map to same space", g);
  984.     }
  985.     holdmat[i] = (temp.m)[0]; 
  986.   }
  987.  
  988. // Build an orthogonal basis using the first two distinct points.  If
  989. // we encounter a duplication, just try the next point.  If that also
  990. // doesn't work, we have an error:
  991.  
  992.   Boolean havedup = FALSE;
  993.   Matrix orthbas(2, numcol);
  994.   Scalar mag = sqrt((holdmat[0] * Transpose(holdmat[0]))[0][0]);
  995.   orthbas[0] = (holdmat[0] / mag)[0];
  996.  
  997.   Matrix candidate = holdmat[1];
  998.   for (int trycount = 1; trycount <= 3; trycount++) {
  999.     if (trycount == 3) { 
  1000.       errh.ErrorExit(sig, "Three points are not distinct", g);
  1001.     }
  1002.     Scalar fac = (candidate * Transpose(orthbas[0]))[0][0];
  1003.     orthbas[1] = (candidate - (fac * orthbas[0]))[0];
  1004.     mag = sqrt((orthbas[1] * Transpose(orthbas[1]))[0][0]);
  1005.     if (fabs(mag) < EPSILON) {
  1006.       candidate = holdmat[2];
  1007.       havedup = TRUE;
  1008.     } else {
  1009.       orthbas[1] = (orthbas[1] / mag)[0];
  1010.       break;
  1011.     }
  1012.   }
  1013.   
  1014. // Project the points into the line specified by the basis we
  1015. // have constructed, and make sure this projection is essentially 
  1016. // superfluous:
  1017.  
  1018.   Matrix ppts(4, 2);
  1019.  
  1020.   for (i = 0; i < 4; i++) {
  1021.     ppts[i] = (holdmat[i] * Transpose(orthbas))[0];
  1022.     Matrix testm = (ppts[i] * orthbas) - holdmat[i];
  1023.     for (int j = 0; j < numcol; j++) {
  1024.       if (fabs(testm[0][j]) > EPSILON) {
  1025.         errh.ErrorExit(sig, "Points are not collinear", g);
  1026.       }
  1027.     }
  1028.   }
  1029.  
  1030. // Use this 2-dim representation to calculate the cross product,
  1031. // using the formulation in Penna & Patterson pg. 163:
  1032.  
  1033.   Scalar fac[4];
  1034.  
  1035.   fac[0] = (ppts[0][0] * ppts[2][1]) - (ppts[2][0] * ppts[0][1]);
  1036.   fac[1] = (ppts[1][0] * ppts[3][1]) - (ppts[3][0] * ppts[1][1]);
  1037.   fac[2] = (ppts[1][0] * ppts[2][1]) - (ppts[2][0] * ppts[1][1]);
  1038.   fac[3] = (ppts[0][0] * ppts[3][1]) - (ppts[3][0] * ppts[0][1]);
  1039.  
  1040. // Catch cases where we do not have enough distinct points:
  1041.  
  1042.   Boolean iszero[4];
  1043.   Boolean haveazero = FALSE;
  1044.   Boolean multizero = FALSE;
  1045.   Boolean allequal = TRUE;
  1046.   for (i = 0; i < 4; i++) {
  1047.     iszero[i] = (fabs(fac[i]) < EPSILON);
  1048.     multizero = (multizero || (haveazero && iszero[i]));
  1049.     haveazero = (haveazero || iszero[i]);
  1050.     allequal = allequal && (fabs(fac[i] - fac[0]) < EPSILON);
  1051.   }
  1052.  
  1053.   if (multizero || allequal) {
  1054.     errh.ErrorExit(sig, "Must have 3 distinct points", g);
  1055.   }
  1056.  
  1057.   AugScalar retval;
  1058.   Scalar denom = fac[2] * fac[3];
  1059.   if (fabs(denom) < EPSILON) {
  1060.     retval = AugScalar(INFINITY);
  1061.   } else {
  1062.     retval = AugScalar((fac[0] * fac[1]) / denom);
  1063.   }
  1064.   return retval;
  1065. }
  1066.  
  1067. // ***********************************************************************
  1068. //
  1069. // Cast a GeOb that is in the one-dimensional space "Reals" into a
  1070. // Scalar.  Cannot use operator Scalar() because resulting automatic
  1071. // casting causes ambiguity in algebraic operations.
  1072. //
  1073.  
  1074. Scalar GeOb::ToScalar(void)
  1075. {
  1076.   if (s == Reals) {
  1077.     return (m[0][0]);
  1078.   } else {
  1079.     errh.ErrorExit("Scalar GeOb::ToScalar(void)",
  1080.                 "Attempted to convert a GeOb not in space Reals to a scalar",
  1081.                 *this);
  1082.   }
  1083. }
  1084.  
  1085. // ***********************************************************************
  1086. // ***********************************************************************
  1087. //
  1088. // Vector Class
  1089. //
  1090. // ***********************************************************************
  1091. // ***********************************************************************
  1092. //
  1093. // Public member functions
  1094. //
  1095. // ***********************************************************************
  1096. //
  1097. // Creates a vector tuple with two elements.
  1098. //
  1099.  
  1100. Vector::Vector(VSpace& ins, GeOb& v1, GeOb& v2)
  1101. {
  1102.   type = VECTOR;
  1103.   *this = Vector(ins, GeObList(v1, v2));
  1104. }
  1105.  
  1106. // ***********************************************************************
  1107. //
  1108. // Creates a vector tuple with three elements.
  1109. //
  1110.  
  1111. Vector::Vector(VSpace& ins, GeOb& v1, GeOb& v2, GeOb& v3)
  1112. {
  1113.   type = VECTOR;
  1114.   *this = Vector(ins, GeObList(v1, v2, v3));
  1115. }
  1116.  
  1117. // ***********************************************************************
  1118. //
  1119. // Need to do memberwise assignment, since matrix class members need
  1120. // to do memberwise assignment:
  1121. //
  1122.  
  1123. Vector& Vector::operator=(Vector& g)
  1124. {
  1125.   holds = g.holds;
  1126.   s = g.s;
  1127.   m = g.m;
  1128.   return (*this);
  1129. }
  1130.  
  1131. // ***********************************************************************
  1132. //
  1133. // Output for debugging:
  1134. //
  1135.  
  1136. void Vector::debug_out(ostream &c, int indent)
  1137. {
  1138.   static char* name = "Vector Object\n";
  1139.   this->data_out(c, indent, name);
  1140.   return;
  1141. }
  1142.  
  1143. // ***********************************************************************
  1144. //
  1145. // Vector creation by specifying coordinates:
  1146. //
  1147.  
  1148. Vector::Vector(VBasis& b, ScalarList& a)
  1149. {
  1150.   static char* sig = "Vector::Vector(VBasis&, ScalarList&)";
  1151.  
  1152.   type = VECTOR;
  1153.   holds = VECTOR;
  1154.   s = b.SpaceOf();
  1155.   
  1156.   if (s.ThisSpaceIs() == TANGENT) { 
  1157.     errh.ErrorExit(sig, "Basis is not in a non-tangent vector space", b, a);
  1158.   }
  1159.  
  1160.   if (a.Length() != s.Dim()) {
  1161.     errh.ErrorExit(sig, "Incorrect number of coordinates specified", b, a);
  1162.   }
  1163.  
  1164.   Matrix hold(1, s.Dim());
  1165.   for (int i = 0; i < s.Dim(); i++) {
  1166.     hold[0][i] = a[i];
  1167.   }
  1168.   m = hold * b.Tostdb();
  1169. }
  1170.  
  1171. // ***********************************************************************
  1172. //
  1173. // Creates a vector tuple in the specified cartesian product space using
  1174. // the objects in the list:
  1175. //
  1176.  
  1177. Vector::Vector(VSpace& ins, GeObList &g)
  1178. {
  1179.   static char* sig = "Vector::Vector(VSpace&, GeObList&)";
  1180.  
  1181.   type = VECTOR;
  1182.   holds = VECTOR;
  1183.   s = ins;
  1184.   m = Matrix(1, ins.MatrixSize());
  1185.  
  1186. // The space must be a non-Tangent Vector Space.  The length of the list
  1187. // must be correct:
  1188.  
  1189.   if (ins.NativeType() != VECTOR) { 
  1190.     errh.ErrorExit(sig, "Space is a tangent vector space", ins);
  1191.   }
  1192.  
  1193.   int n = g.Length();
  1194.   if (n != ins.CPSpaceSize()) {
  1195.     errh.ErrorExit(sig,
  1196.            "Length of list does not match number of elements in space",
  1197.            ins, g);
  1198.   }
  1199.  
  1200. // Go through the list of objects, casting them into vectors or
  1201. // affine vectors as needed.  Report an error if the object cannot
  1202. // be mapped into the required space:
  1203.  
  1204.   for (int j = 0; j < n; j++) {
  1205.  
  1206. //  Get the nth space in the tuple and its coordinate starting slot.
  1207. //  Map the object to the necessary type:
  1208.  
  1209.     Space slotspace = ins[j];
  1210.     int slot = ins.StartSlot(j);
  1211.     GeOb insert = (g[j]).MapTo(slotspace.NativeType());
  1212.  
  1213. //  Check that the spaces match, and if they do, replace the
  1214. //  coordinates in the required slots and go to the next object:
  1215.  
  1216.     if (slotspace == insert.SpaceOf()) {
  1217.       for (int i = 0; i < slotspace.Dim(); i++) {
  1218.         m[0][slot++] = (insert.MatrixRep())[0][i];
  1219.       }
  1220.     } else {
  1221.       errh.ErrorExit(sig, "Mapped GeOb is not in required space",
  1222.                      slotspace, g[j], insert);
  1223.     }
  1224.   }
  1225. }
  1226.  
  1227. // ***********************************************************************
  1228.  
  1229. Boolean Vector::CanMapTo(GeObType t)
  1230. {
  1231.   Boolean rv;
  1232.  
  1233. // Quick kill:
  1234.  
  1235.   if (t == holds) return (TRUE);
  1236.  
  1237.   if ((s.ThisSpaceIs() == TANG_DUAL) || (s.ThisSpaceIs() == LIN_DUAL)) {
  1238.     rv = ((t == VEC_EC) && !(this->IsZeroVector()));
  1239.   } else {
  1240.     switch (t) {
  1241.       case AFF_POINT:
  1242.         rv = (s.StdAffineSubset()).IsIn(*this);
  1243.         break;
  1244.       case AFF_VECTOR:
  1245.     rv = ((s.LinearMapToRef(TANGENT)).Domain()).IsIn(*this);
  1246.         break;
  1247.       case AFF_VEC_EC:
  1248.     rv = ((s.LinearMapToRef(TANGENT)).Domain()).IsIn(*this) &&
  1249.               !(this->IsZeroVector());;
  1250.         break;
  1251.       case VEC_EC:
  1252.     rv = !(this->IsZeroVector());
  1253.         break;
  1254.       case PROJ_POINT:
  1255.     rv = (s.StdAffineSubset()).IsIn(*this);
  1256.         break;
  1257.       case VECTOR:
  1258.       case NULL_GEOB:
  1259.       case ANY_GEOB:
  1260.       default:
  1261.         errh.ErrorExit("Boolean Vector::CanMapTo(GeObType)",
  1262.                "Illegal GeOb Type",
  1263.                    ErrType("Type specified =", t, GEOBTYPES));
  1264.         break;
  1265.     }
  1266.   }
  1267.   return (rv);
  1268. }
  1269.  
  1270. // ***********************************************************************
  1271.  
  1272. GeOb Vector::MapTo(GeObType t)
  1273. {
  1274.   static char* sig = "GeOb Vector::MapTo(GeObType)";
  1275.   GeOb rv;
  1276.  
  1277. // Quick kill:
  1278.  
  1279.   if (t == holds) return (*this);
  1280.  
  1281.   if ((s.ThisSpaceIs() == TANG_DUAL) || (s.ThisSpaceIs() == LIN_DUAL)) {
  1282.     if ((t == VEC_EC) && !(this->IsZeroVector())) {
  1283.       rv = GeOb(VEC_EC, s, m);
  1284.     } else {
  1285.       errh.ErrorExit(sig, "Impossible Map request", *this,
  1286.                      ErrType("Mapping from:", holds, GEOBTYPES),
  1287.                      ErrType("Mapping to:", t, GEOBTYPES));
  1288.     }
  1289.   } else {
  1290.     switch (t) {
  1291.       case AFF_POINT:
  1292.         rv = (s.AffineMapToRef(AFFINE))(*this);
  1293.         break;
  1294.       case AFF_VECTOR:
  1295.         rv = (s.LinearMapToRef(TANGENT))(*this);
  1296.         break;
  1297.       case AFF_VEC_EC:
  1298.         if (!this->IsZeroVector()) {
  1299.           rv = (s.LinearMapToRef(TANGENT))(*this);
  1300.           rv = GeOb(AFF_VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  1301.         } else {
  1302.           errh.ErrorExit(sig, "Can't map zero vector to affine vec. eq. class",
  1303.              *this);
  1304.         }
  1305.         break;
  1306.       case VEC_EC:
  1307.         if (!this->IsZeroVector()) {
  1308.           rv = GeOb(VEC_EC, s, m);
  1309.         } else {
  1310.           errh.ErrorExit(sig, "Can't map zero vector to vec. eq. class",
  1311.              *this);
  1312.         }
  1313.         break;
  1314.       case PROJ_POINT:
  1315.         rv = (s.AffineMapToRef(PROJECT_COMP))(*this);
  1316.         break;
  1317.       case VECTOR:
  1318.       case NULL_GEOB:
  1319.       case ANY_GEOB:
  1320.       default:
  1321.         errh.ErrorExit(sig, "Illegal GeOb type",
  1322.                     ErrType("Type specified =", t, GEOBTYPES));
  1323.         break;
  1324.     }
  1325.   }
  1326.   return (rv);
  1327. }
  1328.  
  1329. // ***********************************************************************
  1330. // ***********************************************************************
  1331. //
  1332. // AVector Class
  1333. //
  1334. // ***********************************************************************
  1335. // ***********************************************************************
  1336. //
  1337. // Public member functions
  1338. //
  1339. // ***********************************************************************
  1340.  
  1341. AVector::AVector(VSpace& ins, GeOb& v1, GeOb& v2)
  1342. {
  1343.   type = AFF_VECTOR;
  1344.   *this = AVector(ins, GeObList(v1, v2));
  1345. }
  1346.  
  1347. // ***********************************************************************
  1348.  
  1349. AVector::AVector(VSpace& ins, GeOb& v1, GeOb& v2, GeOb& v3)
  1350. {
  1351.   type = AFF_VECTOR;
  1352.   *this = AVector(ins, GeObList(v1, v2, v3));
  1353. }
  1354.  
  1355. // ***********************************************************************
  1356. //
  1357. // Need to do memberwise assignment, since matrix class members need
  1358. // to do memberwise assignment:
  1359. //
  1360.  
  1361. AVector& AVector::operator=(AVector& g)
  1362. {
  1363.   holds = g.holds;
  1364.   s = g.s;
  1365.   m = g.m;
  1366.   return (*this);
  1367. }
  1368.  
  1369. // ***********************************************************************
  1370. //
  1371. // Output for debugging:
  1372. //
  1373.  
  1374. void AVector::debug_out(ostream &c, int indent)
  1375. {
  1376.   static char* name = "AVector Object\n";
  1377.   this->data_out(c, indent, name);
  1378.   return;
  1379. }
  1380.  
  1381. // ***********************************************************************
  1382. //
  1383. // Vector creation by specifying coordinates wrt some basis:
  1384. //
  1385.  
  1386. AVector::AVector(Basis& b, ScalarList& a)
  1387. {
  1388.   static char* sig = "AVector::AVector(Basis&, ScalarList&)";
  1389.   int n;
  1390.   Matrix dropdown;
  1391.  
  1392.   type = AFF_VECTOR;
  1393.   holds = AFF_VECTOR;
  1394.  
  1395.   switch (b.Holds()) {
  1396.     case FRAME:
  1397.     case SIMPLEX:
  1398.       s = (b.SpaceOf()).GetSpace(TANGENT);
  1399.       n = b.SpaceOf().MatrixSize();
  1400.       if (a.Length() != n) {
  1401.         errh.ErrorExit(sig, "Incorrect number of coordinates specified", b, a);
  1402.       }
  1403.       dropdown = Transpose((b.SpaceOf()).HoistTangent());
  1404.       if (b.Holds() == SIMPLEX) {
  1405.         Scalar sum = 0.0;
  1406.         for (int i = 0; i < n; i++) {
  1407.           sum += a[i];
  1408.         }
  1409.         if (fabs(sum) > EPSILON) {
  1410.           errh.ErrorExit(sig, "Scalars must sum to 0.0 for a simplex",
  1411.                          ErrVal("Sum = ", sum), a, b);
  1412.         }
  1413.       } else if (fabs(a[n - 1]) > EPSILON) {
  1414.         errh.ErrorExit(sig, "Last scalar must equal 0.0 for a frame",
  1415.                        ErrVal("Value = ", a[n - 1]), b);
  1416.       }
  1417.       break;
  1418.     case VBASIS:
  1419.       s = b.SpaceOf();
  1420.       n = s.MatrixSize();
  1421.       if (a.Length() != n) {
  1422.         errh.ErrorExit(sig, "Incorrect number of coordinates specified", b, a);
  1423.       }
  1424.       dropdown = IdentityMatrix(n);
  1425.       break;
  1426.     case HFRAME:
  1427.     case NULL_BASIS:
  1428.     case ANY_BASIS:
  1429.       default:
  1430.       errh.ErrorExit(sig, "Illegal basis type", b);
  1431.       break;
  1432.   }
  1433.  
  1434. // If a Frame or Simplex is used, we need to change the matrix representation
  1435. // down to the tangent space:
  1436.  
  1437.   Matrix hold(1, n);
  1438.   for (int i = 0; i < n; i++) {
  1439.     hold[0][i] = a[i];
  1440.   }
  1441.   m = hold * b.Tostdb() * dropdown;
  1442. }
  1443.  
  1444. // ***********************************************************************
  1445.  
  1446. AVector::AVector(VSpace& ins, GeObList& g)
  1447. {
  1448.   static char* sig = "AVector::AVector(VSpace&, GeObList&)";
  1449.  
  1450.   type = AFF_VECTOR;
  1451.   holds = AFF_VECTOR;
  1452.   s = ins;
  1453.   m = Matrix(1, ins.MatrixSize());
  1454.  
  1455. // The space must be a Tangent Vector Space.  The length of the list
  1456. // must be correct:
  1457.  
  1458.   if (ins.ThisSpaceIs() != TANGENT) { 
  1459.     errh.ErrorExit(sig, "Space is not a tangent vector space", ins);
  1460.   }
  1461.  
  1462.   int n = g.Length();
  1463.   if (n != ins.CPSpaceSize()) {
  1464.     errh.ErrorExit(sig,
  1465.            "Length of list does not match number of elements in space",
  1466.                    ins, g);
  1467.   }
  1468.  
  1469. // Go through the list of objects, casting them into affine vectors as
  1470. // needed.  Report an error if the object cannot be mapped into the 
  1471. // required space:
  1472.  
  1473.   for (int j = 0; j < n; j++) {
  1474.  
  1475. //  Get the nth space in the tuple and its coordinate starting slot.
  1476. //  Map the object to the necessary type:
  1477.  
  1478.     Space slotspace = ins[j];
  1479.     int slot = ins.StartSlot(j);
  1480.     GeOb insert = (g[j]).MapTo(slotspace.NativeType());
  1481.  
  1482. //  Check that the spaces match, and if they do, replace the
  1483. //  coordinates in the required slots and go to the next object:
  1484.  
  1485.     if (slotspace == insert.SpaceOf()) {
  1486.       for (int i = 0; i < slotspace.Dim(); i++) {
  1487.         m[0][slot++] = (insert.MatrixRep())[0][i];
  1488.       }
  1489.     } else {
  1490.       errh.ErrorExit(sig, "Mapped GeOb is not in required space",
  1491.                       slotspace, g[j], insert);
  1492.     }
  1493.   }
  1494. }
  1495.  
  1496. // ***********************************************************************
  1497. //
  1498. // Returns TRUE iff this AVector can be cast to the specified object:
  1499. //
  1500.  
  1501. Boolean AVector::CanMapTo(GeObType t)
  1502. {
  1503.   Boolean rv;
  1504.  
  1505. // Quick kill:
  1506.  
  1507.   if (t == holds) return (TRUE);
  1508.  
  1509.   switch (t) {
  1510.     case VECTOR:
  1511.       rv = TRUE;
  1512.       break;
  1513.     case AFF_VEC_EC:
  1514.     case VEC_EC:
  1515.       rv = !(this->IsZeroVector());
  1516.       break;
  1517.     case AFF_POINT:
  1518.     case PROJ_POINT:
  1519.       rv = FALSE;
  1520.       break;
  1521.     case AFF_VECTOR:
  1522.     case NULL_GEOB:
  1523.     case ANY_GEOB:
  1524.     default:
  1525.       errh.ErrorExit("Boolean AVector::CanMapTo(GeObType)", "Illegal GeOb Type",
  1526.                  ErrType("Type specified =", t, GEOBTYPES));
  1527.       break;
  1528.   }
  1529.   return (rv);
  1530. }
  1531.  
  1532. // ***********************************************************************
  1533. //
  1534. // Map the AVector to another type:
  1535. //
  1536.  
  1537. GeOb AVector::MapTo(GeObType t)
  1538. {
  1539.   static char* sig = "GeOb AVector::MapTo(GeObType)";
  1540.   GeOb rv;
  1541.  
  1542. // Quick kill:
  1543.  
  1544.   if (t == holds) return (*this);
  1545.  
  1546.   switch (t) {
  1547.     case VECTOR:
  1548.       rv = (s.LinearMapToRef(LINEARIZATION))(*this);
  1549.       break;
  1550.     case VEC_EC:
  1551.       if (!this->IsZeroVector()) {
  1552.         rv = (s.LinearMapToRef(LINEARIZATION))(*this);
  1553.         rv = GeOb(VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  1554.       } else {
  1555.         errh.ErrorExit(sig, "Can't map zero vector to vector eq. class",
  1556.                *this);
  1557.       }
  1558.       break;
  1559.     case AFF_VEC_EC:
  1560.       if (!this->IsZeroVector()) {
  1561.         rv = GeOb(AFF_VEC_EC, s, m);
  1562.       } else {
  1563.         errh.ErrorExit(sig, "Can't map zero vector to affine vec. eq. class",
  1564.                *this);
  1565.       } 
  1566.       break;
  1567.     case AFF_POINT:
  1568.     case PROJ_POINT:
  1569.       errh.ErrorExit(sig, "Impossible map request", *this,
  1570.                  ErrType("Mapping from:", holds, GEOBTYPES),
  1571.                  ErrType("Mapping to:", t, GEOBTYPES));
  1572.       break;
  1573.     case AFF_VECTOR:
  1574.     case NULL_GEOB:
  1575.     case ANY_GEOB:
  1576.     default:
  1577.       errh.ErrorExit(sig, "Illegal GeOb type",
  1578.                  ErrType("Type specified =", t, GEOBTYPES));
  1579.       break;
  1580.   }
  1581.   return (rv);
  1582. }
  1583.  
  1584. // ***********************************************************************
  1585. // ***********************************************************************
  1586. //
  1587. // VectorEC Class
  1588. //
  1589. // ***********************************************************************
  1590. //
  1591. // Public member functions
  1592. //
  1593. // ***********************************************************************
  1594. //
  1595. // Need to do memberwise assignment, since Matrix class members need
  1596. // to do memberwise assignment:
  1597. //
  1598.  
  1599. VectorEC& VectorEC::operator=(VectorEC& g)
  1600. {
  1601.   holds = g.holds;
  1602.   s = g.s;
  1603.   m = g.m;
  1604.   return (*this);
  1605. }
  1606.  
  1607. // ***********************************************************************
  1608. //
  1609. // Output for debugging:
  1610. //
  1611.  
  1612. void VectorEC::debug_out(ostream &c, int indent)
  1613. {
  1614.   static char* name = "VectorEC Object\n";
  1615.   this->data_out(c, indent, name);
  1616.   return;
  1617. }
  1618.  
  1619. // ***********************************************************************
  1620. //
  1621. // Returns TRUE iff this VectorEC can be cast to the specified object:
  1622. //
  1623.  
  1624. Boolean VectorEC::CanMapTo(GeObType t)
  1625. {
  1626.   Boolean rv;
  1627.  
  1628. // Quick kill:
  1629.  
  1630.   if (t == holds) return (TRUE);
  1631.  
  1632.   if ((s.ThisSpaceIs() == TANG_DUAL) || (s.ThisSpaceIs() == LIN_DUAL)) {
  1633.     rv = (t == VECTOR);
  1634.   } else {
  1635.     switch (t) {
  1636.       case AFF_POINT:
  1637.       case AFF_VECTOR:
  1638.       case AFF_VEC_EC:
  1639.         {
  1640.         GeOb test = this->MapTo(VECTOR);
  1641.         rv = ((s.LinearMapToRef(TANGENT)).Domain()).IsIn(test);
  1642.         if (t == AFF_POINT) rv = !rv;
  1643.     } 
  1644.         break;
  1645.       case VECTOR:
  1646.       case PROJ_POINT:
  1647.         rv = TRUE;
  1648.         break;
  1649.       case VEC_EC:
  1650.       case NULL_GEOB:
  1651.       case ANY_GEOB:
  1652.       default:
  1653.         errh.ErrorExit("Boolean VectorEC::CanMapTo(GeObType)",
  1654.                "Illegal GeOb type",
  1655.                    ErrType("Type specified =", t, GEOBTYPES));
  1656.         break;
  1657.     }
  1658.   }
  1659.   return (rv);
  1660. }
  1661.  
  1662. // ***********************************************************************
  1663. //
  1664. // Map this VectorEC to another object:
  1665. //
  1666.  
  1667. GeOb VectorEC::MapTo(GeObType t)
  1668. {
  1669.   static char* sig = "GeOb VectorEC::MapTo(GeObType)";
  1670.   GeOb rv;
  1671.   Scalar coeff;
  1672.  
  1673. // Quick kill:
  1674.  
  1675.   if (t == holds) return (*this);
  1676.  
  1677.   if ((s.ThisSpaceIs() == TANG_DUAL) || (s.ThisSpaceIs() == LIN_DUAL)) {
  1678.     if (t == VECTOR) {
  1679.       coeff = randval();
  1680.       rv = GeOb(VECTOR, s, m * coeff);
  1681.     } else {
  1682.       errh.ErrorExit(sig, "Impossible map request", *this,
  1683.                  ErrType("Mapping from:", holds, GEOBTYPES),
  1684.                  ErrType("Mapping to:", t, GEOBTYPES));
  1685.     }
  1686.   } else {
  1687.     switch (t) {
  1688.       case AFF_POINT:
  1689. //      Go through the back door, by mapping into projective space first!
  1690.         rv = (s.ProjectiveMapToRef(PROJECT_COMP))(*this);
  1691.         rv = ((rv.SpaceOf()).AffineMapToRef(AFFINE))(rv);
  1692.         break;
  1693.       case AFF_VECTOR:
  1694.         coeff = randval();
  1695.         rv = GeOb(VECTOR, s, m * coeff);
  1696.         rv = (s.LinearMapToRef(TANGENT))(rv);
  1697.         break;
  1698.       case AFF_VEC_EC:
  1699.         rv = GeOb(VECTOR, s, m);
  1700.         rv = (s.LinearMapToRef(TANGENT))(rv);
  1701.         rv = GeOb(AFF_VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  1702.         break;
  1703.       case VECTOR:
  1704.         coeff = randval();
  1705.         rv = GeOb(VECTOR, s, m * coeff);
  1706.         break;
  1707.       case PROJ_POINT:
  1708.         rv = (s.ProjectiveMapToRef(PROJECT_COMP))(*this);
  1709.         break;
  1710.       case VEC_EC:
  1711.       case NULL_GEOB:
  1712.       case ANY_GEOB:
  1713.       default:
  1714.         errh.ErrorExit(sig, "Illegal GeOb type",
  1715.                        ErrType("Type specified =", t, GEOBTYPES));
  1716.         break;
  1717.     }
  1718.   }
  1719.   return (rv);
  1720. }
  1721.  
  1722. // ***********************************************************************
  1723. // ***********************************************************************
  1724. //
  1725. // AVectorEC Class
  1726. //
  1727. // ***********************************************************************
  1728. // ***********************************************************************
  1729. //
  1730. // Public member functions
  1731. //
  1732. // ***********************************************************************
  1733. //
  1734. // Need to do memberwise assignment, since Matrix class members need
  1735. // to do memberwise assignment:
  1736. //
  1737.  
  1738. AVectorEC& AVectorEC::operator=(AVectorEC& g)
  1739. {
  1740.   holds = g.holds;
  1741.   s = g.s;
  1742.   m = g.m;
  1743.   return (*this);
  1744. }
  1745.  
  1746. // ***********************************************************************
  1747. //
  1748. // Output for debugging:
  1749. //
  1750.  
  1751. void AVectorEC::debug_out(ostream &c, int indent)
  1752. {
  1753.   static char* name = "AVectorEC Object\n";
  1754.   this->data_out(c, indent, name);
  1755.   return;
  1756. }
  1757.  
  1758. // ***********************************************************************
  1759. //
  1760. // Returns TRUE iff this AVectorEC can be cast to the specified object:
  1761. //
  1762.  
  1763. Boolean AVectorEC::CanMapTo(GeObType t)
  1764. {
  1765.   Boolean rv;
  1766.  
  1767. // Quick kill:
  1768.  
  1769.   if (t == holds) return (TRUE);
  1770.  
  1771.   switch (t) {
  1772.     case AFF_VECTOR:
  1773.     case VEC_EC:
  1774.     case VECTOR:
  1775.     case PROJ_POINT:
  1776.       rv = TRUE;
  1777.       break;
  1778.     case AFF_POINT:
  1779.       rv = FALSE;
  1780.       break;
  1781.     case AFF_VEC_EC:
  1782.     case NULL_GEOB:
  1783.     case ANY_GEOB:
  1784.     default:
  1785.       errh.ErrorExit("Boolean AVectorEC::CanMapTo(GeObType)",
  1786.              "Illegal GeOb type",
  1787.                  ErrType("Type specified =", t, GEOBTYPES));
  1788.       break;
  1789.   }
  1790.   return (rv);
  1791. }
  1792.  
  1793. // ***********************************************************************
  1794. //
  1795. // This is the way to get the GeOb that results from standard mappings:
  1796. //
  1797.  
  1798. GeOb AVectorEC::MapTo(GeObType t)
  1799. {
  1800.   static char* sig = "GeOb AVectorEC::MapTo(GeObType)";
  1801.   GeOb rv;
  1802.   Scalar coeff;
  1803.  
  1804. // Quick kill:
  1805.  
  1806.   if (t == holds) return (*this);
  1807.  
  1808.   switch (t) {
  1809.     case VEC_EC:
  1810.       rv = GeOb(AFF_VECTOR, s, m);
  1811.       rv = (s.LinearMapToRef(LINEARIZATION))(rv);
  1812.       rv = GeOb(VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  1813.       break;
  1814.     case VECTOR:
  1815.       coeff = randval();
  1816.       rv = GeOb(AFF_VECTOR, s, m * coeff);
  1817.       rv = (s.LinearMapToRef(LINEARIZATION))(rv);
  1818.       break;
  1819.     case AFF_VECTOR:
  1820.       coeff = randval();
  1821.       rv = GeOb(AFF_VECTOR, s, m * coeff);
  1822.       break;
  1823.     case PROJ_POINT:
  1824.       rv = GeOb(AFF_VECTOR, s, m);
  1825.       rv = (s.LinearMapToRef(LINEARIZATION))(rv);
  1826.       rv = GeOb(VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  1827.       rv = (rv.SpaceOf().ProjectiveMapToRef(PROJECT_COMP))(rv);
  1828.       break;
  1829.     case AFF_POINT:
  1830.       errh.ErrorExit(sig, "Impossible map request", *this,
  1831.                  ErrType("Mapping from:", holds, GEOBTYPES),
  1832.                  ErrType("Mapping to:", t, GEOBTYPES));
  1833.       break;
  1834.     case AFF_VEC_EC:
  1835.     case NULL_GEOB:
  1836.     case ANY_GEOB:
  1837.     default:
  1838.       errh.ErrorExit(sig, "Illegal GeOb type",
  1839.                  ErrType("Type specified =", t, GEOBTYPES));
  1840.       break;
  1841.   }
  1842.   return (rv);
  1843. }
  1844.  
  1845. // ***********************************************************************
  1846. // ***********************************************************************
  1847. //
  1848. // APoint Class
  1849. //
  1850. // ***********************************************************************
  1851. // ***********************************************************************
  1852. //
  1853. // Public member functions
  1854. //
  1855. // ***********************************************************************
  1856.  
  1857. APoint::APoint(ASpace& ins, GeOb& v1, GeOb& v2)
  1858. {
  1859.   type = AFF_POINT;
  1860.   *this = APoint(ins, GeObList(v1, v2));
  1861. }
  1862.  
  1863. // ***********************************************************************
  1864.  
  1865. APoint::APoint(ASpace &ins, GeOb &v1, GeOb &v2, GeOb &v3)
  1866. {
  1867.   type = AFF_POINT;
  1868.   *this = APoint(ins, GeObList(v1, v2, v3));
  1869. }
  1870.  
  1871. // ***********************************************************************
  1872. //
  1873. // Need to do memberwise assignment, since Matrix class members need
  1874. // to do memberwise assignment:
  1875. //
  1876.  
  1877. APoint& APoint::operator=(APoint& g)
  1878. {
  1879.   holds = g.holds;
  1880.   s = g.s;
  1881.   m = g.m;
  1882.   return (*this);
  1883. }
  1884.  
  1885. // ***********************************************************************
  1886. //
  1887. // Output for debugging:
  1888. //
  1889.  
  1890. void APoint::debug_out(ostream &c, int indent)
  1891. {
  1892.   static char* name = "APoint Object\n";
  1893.   this->data_out(c, indent, name);
  1894.   return;
  1895. }
  1896.  
  1897. // ***********************************************************************
  1898. //
  1899. // APoint creation by specifying coordinates wrt some basis:
  1900. //
  1901.  
  1902. APoint::APoint(Basis& b, ScalarList& a)
  1903. {
  1904.   static char* sig = "APoint::APoint(Basis&, ScalarList&)";
  1905.  
  1906.   holds = AFF_POINT;
  1907.   type = AFF_POINT;
  1908.   s = b.SpaceOf();
  1909.   int n = s.MatrixSize();
  1910.  
  1911.   if (a.Length() != n) {
  1912.     errh.ErrorExit(sig, "Incorrect number of coordinates specified", b, a);
  1913.   }
  1914.  
  1915.   if (b.Holds() == FRAME) {
  1916.     if (fabs(1.0 - a[n - 1]) > EPSILON) {
  1917.       errh.ErrorExit(sig, "Last scalar must equal 1.0 for a frame",
  1918.                      ErrVal("Value = ", a[n - 1]), b);
  1919.     }
  1920.   } else if (b.Holds() == SIMPLEX) {
  1921.     Scalar sum = 0.0;
  1922.     for (int i = 0; i < n; i++) {
  1923.       sum += a[i];
  1924.     }
  1925.     if (fabs(1.0 - sum) > EPSILON) {
  1926.       errh.ErrorExit(sig, "Scalars must sum to 1.0 for a simplex",
  1927.                      ErrVal("Sum = ", sum), a, b);
  1928.     }
  1929.   } else {
  1930.     errh.ErrorExit(sig, "Illegal basis type", b);
  1931.   }
  1932.  
  1933.   Matrix hold(1, n);
  1934.   for (int i = 0; i < n; i++) {
  1935.     hold[0][i] = a[i];
  1936.   }
  1937.   m = hold * b.Tostdb();
  1938. }
  1939.  
  1940. // ***********************************************************************
  1941. //
  1942. // Create a Point tuple:
  1943. //
  1944.  
  1945. APoint::APoint(ASpace& ins, GeObList& g)
  1946. {
  1947.   static char* sig = "APoint::APoint(ASpace&, GeObList&)";
  1948.   int slot;
  1949.  
  1950.   type = AFF_POINT;
  1951.   holds = AFF_POINT;
  1952.   s = ins;
  1953.   m = Matrix(1, ins.MatrixSize());
  1954.  
  1955. // The space must be an Affine Space.  The length of the list must
  1956. // be correct:
  1957.  
  1958.   int n = g.Length();
  1959.   if (n != ins.CPSpaceSize()) {
  1960.     errh.ErrorExit(sig,
  1961.            "Length of list does not match number of elements in space",
  1962.            ins, g);
  1963.   }
  1964.  
  1965. // Go through the list of objects, casting them into points as needed.
  1966. // Report an error if the object cannot be mapped into the required space:
  1967.  
  1968.   for (int j = 0; j < n; j++) {
  1969.  
  1970. //  Get the nth space in the tuple and its coordinate starting slot.
  1971. //  Map the object to the necessary type:
  1972.  
  1973.     Space slotspace = ins[j];
  1974.     slot = ins.StartSlot(j);
  1975.     GeOb insert = (g[j]).MapTo(AFF_POINT);
  1976.  
  1977. //  Check that the spaces match, and if they do, replace the
  1978. //  coordinates in the required slots and go to the next object:
  1979.  
  1980.     if (slotspace == insert.SpaceOf()) {
  1981.       for (int i = 0; i < slotspace.Dim(); i++) {
  1982.         m[0][slot++] = (insert.MatrixRep())[0][i];
  1983.       }
  1984.     } else {
  1985.       errh.ErrorExit(sig, "Mapped GeOb is not in required space",
  1986.                      slotspace, g[j], insert);
  1987.     }
  1988.   }
  1989.  
  1990. // The last slot must be 1.0:
  1991.  
  1992.   m[0][slot] = 1.0;
  1993. }
  1994.  
  1995. // ***********************************************************************
  1996. //
  1997. // Returns TRUE iff this APoint can be mapped to the specified type:
  1998. //
  1999.  
  2000. Boolean APoint::CanMapTo(GeObType t)
  2001. {
  2002.   Boolean rv;
  2003.  
  2004. // Quick kill:
  2005.  
  2006.   if (t == holds) return (TRUE);
  2007.  
  2008.   switch (t) {
  2009.     case AFF_VEC_EC:
  2010.     case AFF_VECTOR:
  2011.       rv = FALSE;
  2012.       break;
  2013.     case VECTOR:
  2014.     case VEC_EC:
  2015.     case PROJ_POINT:
  2016.       rv = TRUE;
  2017.       break;
  2018.     case AFF_POINT:
  2019.     case NULL_GEOB:
  2020.     case ANY_GEOB:
  2021.     default:
  2022.       errh.ErrorExit("Boolean APoint::CanMapTo(GeObType)", "Illegal GeOb type",
  2023.                  ErrType("Type specified =", t, GEOBTYPES));
  2024.       break;
  2025.     }
  2026.   return (rv);
  2027. }
  2028.  
  2029. // ***********************************************************************
  2030. //
  2031. // Maps this object to the specified type:
  2032. //
  2033.  
  2034. GeOb APoint::MapTo(GeObType t)
  2035. {
  2036.   static char* sig = "GeOb APoint::MapTo(GeObType)";
  2037.   GeOb rv;
  2038.  
  2039. // Quick kill:
  2040.  
  2041.   if (t == holds) return (*this);
  2042.  
  2043.   switch (t) {
  2044.     case VECTOR:
  2045.       rv = (s.AffineMapToRef(LINEARIZATION))(*this);
  2046.       break;
  2047.     case VEC_EC:
  2048.       rv = (s.AffineMapToRef(LINEARIZATION))(*this);
  2049.       rv = GeOb(VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  2050.       break;
  2051.     case PROJ_POINT:
  2052.       rv = (s.AffineMapToRef(PROJECT_COMP))(*this);
  2053.       break;
  2054.     case AFF_VEC_EC:
  2055.     case AFF_VECTOR:
  2056.       errh.ErrorExit(sig, "Impossible map request", *this,
  2057.                  ErrType("Mapping from:", holds, GEOBTYPES),
  2058.                  ErrType("Mapping to:", t, GEOBTYPES));
  2059.       break;
  2060.     case AFF_POINT:
  2061.     case NULL_GEOB:
  2062.     case ANY_GEOB:
  2063.     default:
  2064.       errh.ErrorExit(sig, "Illegal GeOb type",
  2065.                  ErrType("Type specified =", t, GEOBTYPES));
  2066.       break;
  2067.   }
  2068.   return (rv);
  2069. }
  2070.  
  2071. // ***********************************************************************
  2072. // ***********************************************************************
  2073. //
  2074. // PPoint Class
  2075. // 
  2076. // ***********************************************************************
  2077. // ***********************************************************************
  2078. //
  2079. // Public member functions
  2080. //
  2081. // ***********************************************************************
  2082. //
  2083. // Need to do memberwise assignment, since Matrix class members need
  2084. // to do memberwise assignment:
  2085. //
  2086.  
  2087. PPoint& PPoint::operator=(PPoint& g)
  2088. {
  2089.   holds = g.holds;
  2090.   s = g.s;
  2091.   m = g.m;
  2092.   return (*this);
  2093. }
  2094.  
  2095. // ***********************************************************************
  2096. //
  2097. // Output for debugging:
  2098. //
  2099.  
  2100. void PPoint::debug_out(ostream &c, int indent)
  2101. {
  2102.   static char* name = "PPoint Object\n";
  2103.   this->data_out(c, indent, name);
  2104.   return;
  2105. }
  2106.  
  2107. // ***********************************************************************
  2108. //
  2109. // Projective Point creation by specifying homogeneous coordinates:
  2110. //
  2111.  
  2112. PPoint::PPoint(HFrame& b, ScalarList& a)
  2113. {
  2114.   type = PROJ_POINT;
  2115.   holds = PROJ_POINT;
  2116.   s = b.SpaceOf();
  2117.   
  2118.   if (a.Length() != s.MatrixSize()) {
  2119.     errh.ErrorExit("PPoint::PPoint(HFrame&, ScalarList&)",
  2120.                    "Incorrect number of coordinates specified", b, a);
  2121.   }
  2122.  
  2123.   Matrix hold(1, s.MatrixSize());
  2124.   for (int i = 0; i < s.MatrixSize(); i++) {
  2125.     hold[0][i] = a[i];
  2126.   }
  2127.   m = hold * b.Tostdb();
  2128. }
  2129.  
  2130. // ***********************************************************************
  2131. //
  2132. // Returns TRUE iff the PPoint can be mapped to the specified object:
  2133. //
  2134.  
  2135. Boolean PPoint::CanMapTo(GeObType t)
  2136. {
  2137.   Boolean rv;
  2138.  
  2139. // Quick kill:
  2140.  
  2141.   if (t == holds) return (TRUE);
  2142.  
  2143.   switch (t) {
  2144.     case VECTOR:
  2145.     case AFF_POINT:
  2146.       rv = (s.StdAffineSubset()).IsIn(*this);
  2147.       break;
  2148.     case AFF_VEC_EC:
  2149.       rv = !((s.StdAffineSubset()).IsIn(*this));
  2150.       break;
  2151.     case AFF_VECTOR:
  2152.       rv = FALSE;
  2153.       break;
  2154.     case VEC_EC:
  2155.       rv = TRUE;
  2156.       break;
  2157.     case PROJ_POINT:
  2158.     case NULL_GEOB:
  2159.     case ANY_GEOB:
  2160.     default:
  2161.       errh.ErrorExit("Boolean PPoint::CanMapTo(GeObType)", "Illegal GeOb type",
  2162.                  ErrType("Type specified =", t, GEOBTYPES));
  2163.       break;
  2164.   }
  2165.   return (rv);
  2166. }
  2167.  
  2168. // ***********************************************************************
  2169. //
  2170. // This is the way to get the GeOb that results from standard mappings:
  2171. //
  2172.  
  2173. GeOb PPoint::MapTo(GeObType t)
  2174. {
  2175.   static char* sig = "GeOb PPoint::MapTo(GeObType)";
  2176.   GeOb rv;
  2177.  
  2178. // Quick kill:
  2179.  
  2180.   if (t == holds) return (*this);
  2181.  
  2182.   switch (t) {
  2183.     case VECTOR:
  2184.       rv = (s.AffineMapToRef(LINEARIZATION))(*this);
  2185.       break;
  2186.     case AFF_POINT:
  2187.       rv = (s.AffineMapToRef(AFFINE))(*this);
  2188.       break;
  2189.     case VEC_EC:
  2190.       rv = (s.ProjectiveMapToRef(LINEARIZATION))(*this);
  2191.       break;
  2192.     case AFF_VEC_EC:
  2193.       rv = (s.ProjectiveMapToRef(LINEARIZATION))(*this);
  2194.       rv = GeOb(VECTOR, rv.SpaceOf(), rv.MatrixRep());
  2195.       rv = (rv.SpaceOf().LinearMapToRef(TANGENT))(rv);
  2196.       rv = GeOb(AFF_VEC_EC, rv.SpaceOf(), rv.MatrixRep());
  2197.       break;
  2198.     case AFF_VECTOR:
  2199.       errh.ErrorExit(sig, "Impossible Map request", *this,
  2200.                      ErrType("Mapping from:", holds, GEOBTYPES),
  2201.                      ErrType("Mapping to:", t, GEOBTYPES));
  2202.       break;
  2203.     case PROJ_POINT:
  2204.     case NULL_GEOB:
  2205.     case ANY_GEOB:
  2206.     default:
  2207.       errh.ErrorExit(sig, "Illegal GeOb Type",
  2208.                  ErrType("Type specified =", t, GEOBTYPES));
  2209.       break;
  2210.   }
  2211.   return (rv);
  2212. }
  2213.  
  2214. // ***********************************************************************
  2215. // ***********************************************************************
  2216. //
  2217. // Create a "standard" Null GeOb that can be used for building lists
  2218. // of geometric object arguments for partial multimap evaluation:
  2219. //
  2220.  
  2221. GeOb NullGeOb;
  2222.  
  2223. // ***********************************************************************
  2224.  
  2225.  
  2226.  
  2227.  
  2228.